The spectral weight of the Hubbard model through cluster perturbation theory 



D. Scnechal, D. Perez and M. Pioro-Ladriere 
Centre de Recherche en Physique du SoUde et Departement de Physique, 
Universite de Sherbrooke, Sherbrooke, Quebec, Canada JIK 2R1. 
(August 1999) 

We calculate the spectral weight of the one- and two-dimensional Hubbard models, by performing 
exact diagonahzations of finite clusters and treating inter-cluster hopping with perturbation theory. 
Even with relatively modest clusters (e.g. 12 sites), the spectra thus obtained give an accurate 
description of the exact results. Thus, spin-charge separation (i.e. an extended spectral weight 
bounded by singularities) is clearly recognized in the one-dimensional Hubbard model, and so is 
extended spectral weight in the two-dimensional Hubbard model. 



One of the central issues in the theory of strongly cor- 
related electrons is the existence or not of well-defined 
quasiparticles. This question is best addressed by study- 
ing the spectral weight (SW) A(k, w), i.e., the probabifity 
distribution for the energy huj of an electron of wavevec- 
tor k added to, or removed from the system. In a Fermi 
liquid, the SW is dominated by a single quasiparticle 
peak centered at u; = £(k), whose width decreases as 
e(k) approaches the Fermi energy. In a Luttinger liquid, 
the SW is distributed between two singularities associ- 
ated respectively to spin and charge excitations (spinous 
and holons) The hole (i.e., electron-removal) part of 
the SW can be measured by angle-resolved photoemission 
spectroscopy (ARPES), a technique which has improved 
steadily in recent years On the theoretical side, the 

SW is obtained from the one-particle Green function: 



A(k,w)= lim -2Imt;(k,w 



(1) 



and the latter may be approximately evaluated by various 
analytical and numerical methods. 

In this letter we explain a new method for calculat- 
ing A{\<i,uj) in Hubbard-type models, based on a combi- 
nation of exact diagonahzations (ED) of finite clusters 
with strong-coupling perturbation theory |^J|], and ap- 
ply it to the one- and two-dimensional Hubbard models. 
Exact diagonahzations based on the Lanczos algorithm 
are commonly used to evaluate ^(k, cj) [0-0. Unfor- 
tunately, computer memory requirements grow exponen- 
tially with system size and restrict the analysis to small 
clusters (e.g. 16 sites for the Hubbard model). The SW 
thus obtained is the sum of a relatively small number of 
poles and its extended character in the thermodynamic 
limit is difficult to assess. The SW may also be evaluated 
by quantum Monte Carlo (QMC) 0-|l|l: lar ger systems 
may thus be studied (e.g. 64 sites) but the maximum 
entropy method (MEM) used for approximate analytic 
continuation tends to produce smooth SW and may miss 
weak features; moreover, computation time increases as 
the temperature is lowered. The new method we pro- 
pose consists in (i) dividing the lattice into identical iV- 
site clusters; (ii) evaluating - by ED - the one-particle 



Green function Ga.b{z) within a cluster {a,b are lattice 
sites and z a complex frequency) with open boundary 
conditions; (iii) treating the inter-cluster hopping in 
perturbation theory and recovering the Green function 
Q{k,u}). Thus, short-distance effects are treated exactly, 
while long-distance propagation is treated at the RPA 
level. 
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FIG. 1. Dividing the lattice into identical four-site clusters 
for the ID and 2D Hubbard models. The in-cluster hopping 
amplitude is t and the inter-cluster hopping is to- 

Step (i) is illustrated on Fig. |^. We denote by t and 
to the hopping amplitudes within and between clusters, 
respectively {t and may be different a priori, but will be 
identical in practice). Clusters of up to = 12 sites have 
been treated, with various aspect ratios in the 2D case. 
Open (free) boundary conditions must be used. Step (ii) 
proceeds according to the usual Lanczos method . The 
cluster Green function Ga.b{z) is defined as 



Ga.t{z) = {^\ca — TTcl\n) + (r!|4^— -c, |r!) (2) 
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z + H 



where jfi) is the ground state obtained by ED, Ca is the 
electron destruction operator at site a (we drop the spin 
index) and H is the Hamiltonian (including chemical po- 
tential). The two terms in Ga.b correspond respectively 
to electron (G° ^) and hole (Gjj f,) propagation, and are 
calculated separately. In the subspace containing one 
additional electron (with respect to the ground state), 
an accurate tridiagonal representation of H is obtained 
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(typically of dimension ranging from 50 to 250). Effi- 
cient routines for inverting tridiagonal matrices are used 
to evaluate G^^(2) at any desired complex value, and 
likewise for G^fj{z). In usual Lanczos calculations, this 
inversion takes the form of a continued-fraction represen- 
tation. 

Step (iii) demands more explanations. Let Crn,a be the 
electron destruction operator on site a of cluster m (a = 
1, . . . , N). The full system is treated as a superlattice of 
clusters, each cluster being made of N ordinary lattice 
sites; we will work in one dimension for simplicity, but 
a suitable generalization to higher-dimensional lattices 
is readily obtained. The complete Hamiltonian of the 
system may be written as H ~ Hq + V: 



m,n 
a.h 



(3) 



where is, say, the Hubbard Hamiltonian of the m}^ 
cluster 

Hi ^-^Y.Y. (4,a,.c™,b,. + H.c.) 

(a, 6) cr 

(4) 

a 

and V is the nearest-neighbor hopping between adjacent 
clusters: 



-to{Sm,n~'lSa,N^b,l + <^m,n+l^aa^6,Af) (5) 



Of interest is the electron Green function Ga,biQ, z), 
where Q is a superlattice wavevector, and a, b are site 
indices within a cluster. The perturbation V being a 
one-body operator, it may be treated in the formalism 
of Refs [^j6|, wherein a systematic perturbation expan- 
sion was constructed for such terms. The lowest-order 
contribution to this expansion has the RPA form: 



G{z) 



1 - V{Q)G{z) 



(6) 



a,b 



where G{z) is the generalization of the "atomic" Green 
function, now a N x N matrix in the space of site indices. 
Likewise, V{Q) is the reciprocal superlattice representa- 
tion of the hopping (^): 



VaM = -to {e'QSa,N6bA+e~''^Sa,l5l,^N) 



(7) 



Relation may be regarded as a cluster generalization 
of the Hubbard-I approximation jl^] . 

The Green function GaAQ) of Eq. (|) is in a mixed 
representation: real space within a cluster and Fourier 
space between clusters. A true Fourier representation in 
terms of the original reciprocal lattice is preferred. Since 
the cluster decomposition breaks translation invariance. 



Q will depend on two continuous momenta k and k' , iden- 
tical modulo a reciprocal superlattice vector: 



N 

X J2 e-'''^''-^^e^"'^/^ga,b{Nk,z) (8) 



If we set t — to, the s — component (fc = k') is the 
RPA approximation to the full Green function: 



a,b 



(9) 



Eqs (^,^) are then used to calculate the SW. 

The approximation (^) turns out to be exact in the 
absence of interactions {U = 0). In that case. Wick's 
theorem applies and the RPA form is the exact resum- 
mation of the perturbation series, V{Q) being the exact 
self-energy. If we set t = tg, Eq. then yields the 
exact Green function for an infinite system at arbitrary 
wavevector, from the exact Green function G of an finite, 
open cluster. When U ^ 0, Expression is no longer 
exact, but strong interactions tend to cause short-range 
correlations that are incorporated with good accuracy in 
modest-size clusters. Thus short-distance effects are well 
served by the ED within a cluster, and long-distance ef- 
fects by the RPA approximation, making our method 
adequate at intermediate coupling. 
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FIG. 2. Above: A{-k/2,uj) for the ID Hubbard model at 
half-filling on a 12-site ring with ordinary ED. Below: the 
same, but with applying RPA with 12-site clusters. The ex- 
tended character of the SW is manifest. Note: the parameter 
77 of Eq. (^) has been set to 0.03 in order to give delta peaks 
a finite width. 
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We have applied the method just described to the ID 
Hubbard model. Fig. |2| compares the result of an or- 
dinary ED on a 12-site ring with the present approach: 
Whereas the extended character of A{k,uj) - here a sig- 
nature of spin-charge separation - is not clear form the 
ED data alone, it is clearly revealed by the RPA spec- 
tra. In fact, the two branches of the SW can already be 
seen with a two-site cluster (not shown), but more and 
more poles appear in between when the cluster size is 
increased: the actual separation of spin and charge exci- 
tations needs a fair cluster size to occur, and propagation 
between clusters at the RPA level requires the holon and 
spinon to recombine. 
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FIG. 3. Spectral weight of the ID Hubbard model at 
half-filling, for U = 4t, calculated from Eqs. and with 
A'' = 12. Below, density plot of the same data. 

Fig. H illustrates the SW of the ID Hubbard model at 
half-filling with U/t — 4. Most noteworthy are: (i) The 
extended character of the SW, with six branches having 
a clear dispersion, even though most of the weight lies 
near the "quasiparticle band" following approximately 
the —2t cos k free-particle dispersion; (ii) the gap opening 
at A; = 7r/2; (iii) the spinon (A) and holon (B) branches, 
characteristic of a Luttinger liquid with a charge gap 
(Branch D is the mirror of the holon branch with op- 
posite frequency) (iv) the weak, higher-frequency 
band (C), absent from low-energy Luttinger liquid pre- 
dictions; (v) the high-frequency tail (E) near the zone 
boundary. Band C, as well as Bands B and D together, 
disperse with period tt instead of 27r, a signature of local 
AF correlations. A comparison with Fig. Ic of Ref. 0] 
- which illustrates the SW in the U ^ oo limit - is re- 
vealing of the changes brought about by a finite value of 



U: in the U ^ oo limit, just the hole part of the SW 
is present, but the same branches can be found, however 
with comparable relative intensity: Branches D and E 
are the mirror images of branches B and A, respectively. 
The finite value of U weakens considerably the intensity 
of branches C, D, and E. It is also interesting to compare 
Fig. I with Fig. 2 of Ref. Q and Fig. 3 of Ref. g, where 
the same parameters were used. In particular, it is clear 
that the Maximum Entropy Method of Ref. lumps 
the spinon and holon peaks near k = into one broad 
peak. 
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4. Spectral weight of the ID Hubbard model at 
and U = At, calculated from Eqs. ^ and ^ with 
12. Below: density plot of the same data. 



Fig. ^ shows the SW of the ID Hubbard model for the 
same value oi U/t, but away from half- filling, at (n) = |. 
The chemical potential fi = 0.64 was inferred from the 
integrated density of states. The Fermi level crosses the 
main band, causing metallic behavior, but again the SW 
is clearly extended, with a clear weakening of the upper 
Hubbard band: there is a significant transfer of SW from 
high to low frequencies jl7|. Again, the spinon (A) and 
holon (B) branches are clearly identified, this time in a 
gapless Luttin ger liquid. This may be compared with 
Fig. 3 of Ref. ||lj], which corresponds to the same U/t 
ratio, but with (n) — |. 

We have also applied our method to the 2D Hubbard 
model, with various cluster shapes (2 x 6, 3 x 4, 4 x 3 
and 6x2). The spectra obtained from these different 
shapes are very similar, and this reinforces our confidence 
in them, even though the linear dimensions of the clus- 
ters are modest. Fig. |^ illustrates the SW of the 2D 
Hubbard model at half-filling for U = 8t, with a 3 x 4 
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cluster. This is to be compared with Fig. 1 of Ref. 
and Fig. 6 of Ref. In contrast to the ID case, the 
SW is much more concentrated around one peak, but 
its extended character is stih undeniable. Indeed, one is 
tempted to draw an analogy with ID spinon and holon 
branches: the momentum scan F — X — M shares fea- 
tures with the [0, tt] scan in the ID case, except that the 
"spinon" (A) is much weaker than the "holon" (B). The 
same can be said of the diagonal scan F — M (from right 
to left on Fig. |^). Likewise, a high-frequency band (C) 
is visible. Most obvious is the gap opening at the Fermi 
surface, constant along the XY line, a feature that would 
certainly be modified by including a diagonal hopping t' , 
and which demonstrates anyway that nearest-neighbor 
hopping alone cannot account for the ARPES data of in- 
sulating cuprates Q (this was already known for the t — J 
model). Note, that whereas Refs. |Q and ^ both resolve 
two peaks near Point M, the present approach suggests 
an extended SW at that point. The expected antiferro- 
magnetic order of the half-filled 2D Hubbard model is 
not seen here, because of finite cluster size. This order 
would imply a folding of the Brillouin zone, with a cor- 
responding symmetry of the SW following the SDW fit 
illustrated on Fig. H. 
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FIG. 5. Spectral weight of the 2D Hubbard model at 
half-filling, for U — 8t, calculated from Eqs. (|^,^) with a 3 x 4 
cluster. Inset: The wavevector scans. Below: density plot 
of the same data (except for the scan from X to Y) to be 
compared directly with Fig. 1 of Ref. |l3| . The dashed curve 
is the best SDW dispersion, with t = 1.36 and gap A — 2.21. 



Some general remarks are in order. Formulas (|6|j9|) 
are but the first order result of a systematic expan- 
sion (See Ref. 1^ for details). It is difficult to assess the 



convergence of this perturbative expansion, since it de- 
pends certainly on the ratio t/U and on cluster size N. 
We expect nonetheless the method to give better results 
at strong coupling, where short-range effects dominate 
and are thus well accounted for by modest clusters. In- 
deed, the effect of antiferromagnetic correlations are al- 
ready seen with two-site clusters. Going to order t§ in 
strong-coupling perturbation is a way of improving the 
results presented here, but appears quite difficult in prac- 
tice, because of the need to compute numerically exact 
two-particle Green functions on a cluster. The spectra 
presented here are all normalized, up to 1 or 2%. The 
general form (^ guarantees that the continued fraction 
form of the SW will have the correct ffist coefficient, thus 
ensuring its normalization. 

In summary, we have shown how strong-coupling per- 
turbation theory can be used to incorporate long-distance 
effects into ED data which already contain short-distance 
effects exactly. This method allows for a clear recogni- 
tion of spin-charge separation in the ID Hubbard model, 
and of extended SW in the 2D Hubbard model. Further 
applications of this method (NNN hopping, three-band 
Hubbard model, etc.) are under way. 
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